Data Mining on CHUSJ dataset
João Miguel Alves
2023-10-21
Fig.1 - Date of 1st lab results
Weekly temporal distribution of first positive COVID-19 lab results throughout the study period.
FirstDate <-
BD_chusj %>%
select(DateOfFirstPositiveLabResult) %>%
group_by(date = DateOfFirstPositiveLabResult) %>%
summarise(freq = n())
FirstDate_daily <- xts::as.xts(FirstDate$freq, order.by = as.Date(FirstDate$date))
FirstDate_weekly <- xts::apply.weekly(FirstDate_daily, sum)
xts::plot.xts(FirstDate_weekly, col = "cadetblue4",
lwd = 3, main = "Date of first positive lab result - weekly",
yaxis.right = FALSE)
Fig.2 - Radar Plots
Radar plots representing medication usage in subgroups: (A) ICU-admitted patients, (B) type of hospital discharge, and (C) type of MRCS. In each radar, the lines represent distinct patient categories, and the distance from each point to the origin indicates the percentage of patients using drugs from a specific ATC group. For instance, a line extending to the 75% marker under the label “ATC.A” suggests that 75% of those patients are being medicated with drugs from the ATC.A group. It’s important to note that comparisons between lines within the same radar provide insights into discrepancies in medication usage among the subgroups; however, specific values and detailed interpretation should be derived from the study’s context.
#' Define function to revert order of variables in radar plot
clockwise_order <- function(ds){
# ds <-ATCs.lv1
ds_rev <- rev(ds)
ds_rev <- ds_rev[, c(ncol(ds_rev), 1:(ncol(ATCs.lv1)-1))]
return(ds_rev)
}
# Define multiframe
par(mfrow = c(1, 3))
legend.x = -1.3
legend.y = -0.9
# Fig.2A - ATCs VS Intensive Care -----------------------------------------
ATCs.lv1 <- BD_chusj[, c(469, 378:390)]
ATCs.lv1 <- cbind(IntensiveCare = BD_chusj$IntensiveCare,
ATCs.lv1)
IC_yes_count <- ATCs.lv1$IntensiveCare[ATCs.lv1$IntensiveCare == 1] %>% length()
IC_no_count <- ATCs.lv1$IntensiveCare[ATCs.lv1$IntensiveCare == 0] %>% length()
# Calculate the percentage of patients in and out of IC for each ATC group
ATCs.lv1 <- rbind(
# Intensive Care == YES
IC_yes = ATCs.lv1 %>%
filter(IntensiveCare == 1) %>%
{sapply(., sum) / IC_yes_count} %>%
round(digits = 3) %>% t() %>% as.data.frame(),
# Intensive Care == NO
IC_no = ATCs.lv1 %>%
filter(IntensiveCare == 0) %>%
{sapply(., sum) / IC_no_count} %>%
round(digits = 3) %>% t() %>% as.data.frame()
)
ATCs.lv1$IntensiveCare <- NULL
# Find max and round it (ceiling) to the nearest .1 multiple
max.axis <- plyr::round_any(max(ATCs.lv1), 0.1, f = ceiling)
# Prepare data for radar chart
ATCs.lv1 <- rbind(max = rep(max.axis, length(ATCs.lv1)),
min = rep(0, length(ATCs.lv1)),
ATCs.lv1)
# Define colors
colors <- c("indianred2", "forestgreen")
colors_border <- colors
colors_in <- scales::alpha(colors, alpha = 0.2)
# Revert order of columns
ATCs.lv1 <- clockwise_order(ATCs.lv1)
radarchart(ATCs.lv1,
axistype = 1,
# custom polygon
pcol = colors_border, # line color
pfcol = colors_in, # fill-in color
plwd = 2, # line width
plty = 1, # line type
# custom the grid
cglcol = "grey20", # net color
cglty = 3, # net line type
axislabcol = "grey50", # net label color
cglwd = 0.8, # net width
caxislabels = scales::percent(seq(0, max.axis, length.out = 5)),
# custom labels
calcex = 0.8,
vlcex = 0.7,
title = "Relative use of medication in ICU-admitted\nand non-ICU-admitted patients",
sub = paste0("ICU-admitted patients: n=", IC_yes_count, "\n",
"Non-ICU-admitted patients: n=", IC_no_count)
)
legend(x = legend.x,
y = legend.y,
legend = rownames(ATCs.lv1[-c(1,2),]),
bty = "n",
pch = 20,
col = colors_border,
text.col = "grey40",
cex = 0.9,
pt.cex = 2)
text(-1.2, 1.2, "(A)", cex = 1.4, pos = 1, col = "black")
# Fig.2B - ATCs vs Hospital Discharge -------------------------------------
ATCs.lv1 <- BD_chusj[, c(469, 378:390)]
ATCs.lv1 <- cbind(Outcome = BD_chusj$Outcome,
ATCs.lv1)
Outcome_Recovered_count <- ATCs.lv1$Outcome[ATCs.lv1$Outcome == 0] %>% na.omit() %>% length()
Outcome_Transferred_count <- ATCs.lv1$Outcome[ATCs.lv1$Outcome == 2] %>% na.omit() %>% length()
Outcome_Deceased_count <- ATCs.lv1$Outcome[ATCs.lv1$Outcome == 1] %>% na.omit() %>% length()
# Calculate the percentage of patients in and out of IC for each ATC group
ATCs.lv1 <- rbind(
Recovered = ATCs.lv1 %>% filter(Outcome == 0) %>% {sapply(., sum) / Outcome_Recovered_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
Transferred = ATCs.lv1 %>% filter(Outcome == 2) %>% {sapply(., sum) / Outcome_Transferred_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
Deceased = ATCs.lv1 %>% filter(Outcome == 1) %>% {sapply(., sum) / Outcome_Deceased_count} %>% round(digits = 3) %>% t() %>% as.data.frame()
)
ATCs.lv1$Outcome <- NULL
# Find max and round it (ceiling) to the nearest .1 multiple
max.axis <- plyr::round_any(max(ATCs.lv1), 0.1, f = ceiling)
# Prepare data for radar chart
ATCs.lv1 <- rbind(max = rep(max.axis, length(ATCs.lv1)),
min = rep(0, length(ATCs.lv1)),
ATCs.lv1)
# Define colors
colors <- c("forestgreen", "dodgerblue3", "firebrick3")
colors_border <- colors
colors_in <- scales::alpha(colors, alpha = 0.1)
# Revert order of columns
ATCs.lv1 <- clockwise_order(ATCs.lv1)
radarchart(ATCs.lv1,
axistype = 1,
# custom polygon
pcol = colors_border, # line color
pfcol = colors_in, # fill-in color
plwd = 2, # line width
plty = 1, # line type
# custom the grid
cglcol = "grey20", # net color
cglty = 3, # net line type
axislabcol = "grey50", # net label color
cglwd = 0.8, # net width
caxislabels = scales::percent(seq(0, max.axis, length.out = 5)),
# custom labels
calcex = 0.8,
vlcex = 0.7,
title = "Relative use of medication in patients\n based on hospital discharge",
sub = paste0("Patients Recovered: n=", Outcome_Recovered_count, "\n",
"Patients Transfered: n=", Outcome_Transferred_count, "\n",
"Patients Deceased: n=", Outcome_Deceased_count)
)
legend(x = legend.x,
y = legend.y,
legend = rownames(ATCs.lv1[-c(1,2),]),
bty = "n",
pch = 20,
col = colors_border,
text.col = "grey40",
cex = 0.9,
pt.cex = 2)
text(-1.2, 1.2, "(B)", cex = 1.4, pos = 1, col = "black")
# Fig.2C - ATCs vs Respiratory Support ------------------------------------
ATCs.lv1 <- BD_chusj[, c(469, 378:390)]
ATCs.lv1 <- cbind(BD_chusj[25:28],
ATCs.lv1)
Support_ECMO_count <- ATCs.lv1$ECMO[ATCs.lv1$ECMO == 1] %>% na.omit() %>% length()
Support_HFO_count <- ATCs.lv1$HFO[ATCs.lv1$HFO == 1] %>% na.omit() %>% length()
Support_NIV_count <- ATCs.lv1$NIV[ATCs.lv1$NIV == 1] %>% na.omit() %>% length()
Support_IMV_count <- ATCs.lv1$IMV[ATCs.lv1$IMV == 1] %>% na.omit() %>% length()
# Calculate the percentage of patients in and out of IC for each ATC group
ATCs.lv1 <- rbind(
ECMO = ATCs.lv1 %>% filter(ECMO == 1) %>% {sapply(., sum) / Support_ECMO_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
HFO = ATCs.lv1 %>% filter(HFO == 1) %>% {sapply(., sum) / Support_HFO_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
NIV = ATCs.lv1 %>% filter(NIV == 1) %>% {sapply(., sum) / Support_NIV_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
IMV = ATCs.lv1 %>% filter(IMV == 1) %>% {sapply(., sum) / Support_IMV_count} %>% round(digits = 3) %>% t() %>% as.data.frame()
)
ATCs.lv1$ECMO <- ATCs.lv1$HFO <- ATCs.lv1$NIV <- ATCs.lv1$IMV <- NULL
# Find max and round it (ceiling) to the nearest .1 multiple
max.axis <- plyr::round_any(max(ATCs.lv1), 0.1, f = ceiling)
# Prepare data for radar chart
ATCs.lv1 <- rbind(max = rep(max.axis, length(ATCs.lv1)),
min = rep(0, length(ATCs.lv1)),
ATCs.lv1)
# Define colors
colors <- c("firebrick2", "dodgerblue2", "forestgreen", "darkorange")
colors_border <- colors
colors_in <- scales::alpha(colors, alpha = 0.1)
# Revert order of columns
ATCs.lv1 <- clockwise_order(ATCs.lv1)
radarchart(ATCs.lv1,
axistype = 1,
# custom polygon
pcol = colors_border, # line color
pfcol = colors_in, # fill-in color
plwd = 2, # line width
plty = 1, # line type
# custom the grid
cglcol = "grey20", # net color
cglty = 3, # net line type
axislabcol = "grey50", # net label color
cglwd = 0.8, # net width
caxislabels = scales::percent(seq(0, max.axis, length.out = 5)),
# custom labels
calcex = 0.8,
vlcex = 0.7,
title = "Relative use of medication in patients based on\n Mechanical Respiratory and Circulatory Support",
sub = paste0("Patients using support:\n",
"ECMO: n=", Support_ECMO_count, "; ",
"HFO: n=", Support_HFO_count, ";\n",
"NIV: n=", Support_NIV_count, "; ",
"IMV: n=", Support_IMV_count
)
)
legend(x = legend.x,
y = legend.y,
legend = rownames(ATCs.lv1[-c(1,2),]),
bty = "n",
pch = 20,
col = colors_border,
text.col = "grey40",
cex = 0.9,
pt.cex = 2)
text(-1.2, 1.2, "(C)", cex = 1.4, pos = 1, col = "black")
# Clear multiframe
par(mfrow = c(1,1))Fig.3 - ATC families, 1st level
Heatmap illustrating co-frequencies of exposure among 1st level ATC drug groups. Each cell represents the relative co-frequency of patients using two drug groups concurrently. For instance, the cell “ATC.A x ATC.C”, with a value of 0.360, indicates that 36% of the total patients (n=2688) are simultaneously prescribed drugs from ATC groups A and C. The diagonal line, like the cell “ATC.C x ATC.C” showing 0.412, can be interpreted as the relative frequency of the use of that specific ATC drug group on its own, signifying that 41.2% of the total patients are prescribed only drugs from ATC group C. It’s important to note that the heatmap is mirrored, meaning data below the diagonal mirror those above it. Therefore, the cell “ATC.A x ATC.C” has the same value as the cell “ATC.C x ATC.A”.
#' 469 + 378-390: ATC codes - first level
ATC_lv1 <- BD_chusj[, c(469, 378:390)]
# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA, ncol = length(ATC_lv1), nrow = length(ATC_lv1)))
colnames(heatmap.matrix) <- rownames(heatmap.matrix) <- names(ATC_lv1)
for (i in 1:length(ATC_lv1)) {
for (j in i:length(ATC_lv1)) {
heatmap.matrix[i,j] <- heatmap.matrix[j,i] <- nrow(ATC_lv1[,c(i,j)][ATC_lv1[,i] == 1 & ATC_lv1[,j] == 1,])
}
}; rm(i, j)
heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value")
# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(ATC_lv1), digits = 3)
# Plot heapmap
ggplot(heatmap.data, aes(x, y, fill = value)) +
geom_tile(color = "grey95",
lwd = 1,
linetype = 1) +
theme(axis.text.x = element_text(angle = 90),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
subtitle = paste0("Expressed as % of total patients (n = ", nrow(ATC_lv1)," patients)")) +
geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
# scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
guides(fill = guide_colourbar(barwidth = 0.5,
barheight = 20,
title = "%"))
Triangular version 1
Excluding diagonal line
# Reshape data to triangular matrix - version_1
heatmap.data_aux <- reshape2::dcast(heatmap.data, x ~ y)
n <- nrow(heatmap.data_aux)
mat <- matrix(NA, n, n)
mat[row(mat) + col(mat) <= n] <- 1
mat <- mat |> as.data.frame() |> rev()
data_vector <- as.vector(as.matrix(mat)) * as.vector(as.matrix(heatmap.data_aux[, 2:(n+1)]))
heatmap.data_diagonal_v1 <- cbind(heatmap.data_aux$x, as.data.frame(matrix(data_vector, nrow = n)))
names(heatmap.data_diagonal_v1) <- names(heatmap.data_aux)
heatmap.data_diagonal_v1 <- reshape2::melt(heatmap.data_diagonal_v1, value.name = "value", id.vars = "x")
names(heatmap.data_diagonal_v1) <- c("x", "y", "value")
heatmap.data_diagonal_v1 <- heatmap.data_diagonal_v1[complete.cases(heatmap.data_diagonal_v1),]
# Plot triangular heapmap v1
ggplot(heatmap.data_diagonal_v1, aes(x, y, fill = value)) +
theme_light() +
geom_tile(color = "grey95",
lwd = 1,
linetype = 1) +
theme(axis.text.x = element_text(angle = 90),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.border=element_blank(),
panel.grid.major=element_line(color='white')) +
ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
subtitle = paste0("Expressed as % of total patients (n = ", nrow(ATC_lv1)," patients)")) +
geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
# scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
guides(fill = guide_colourbar(barwidth = 0.5,
barheight = 20,
title = "%"))
Triangular version 2
Excluding diagonal line
# Reshape data to triangular matrix - version_2
heatmap.data_aux <- reshape2::dcast(heatmap.data, x ~ y)
n <- nrow(heatmap.data_aux)
mat <- matrix(NA, n, n)
mat[row(mat) + col(mat) > n+1] <- 1
mat <- mat |> as.data.frame() |> rev()
data_vector <- as.vector(as.matrix(mat)) * as.vector(as.matrix(heatmap.data_aux[, 2:(n+1)]))
heatmap.data_diagonal_v2 <- cbind(heatmap.data_aux$x, as.data.frame(matrix(data_vector, nrow = n)))
names(heatmap.data_diagonal_v2) <- names(heatmap.data_aux)
heatmap.data_diagonal_v2 <- reshape2::melt(heatmap.data_diagonal_v2, value.name = "value", id.vars = "x")
names(heatmap.data_diagonal_v2) <- c("x", "y", "value")
heatmap.data_diagonal_v2 <- heatmap.data_diagonal_v2[complete.cases(heatmap.data_diagonal_v2),]
# Plot triangular heapmap v2
ggplot(heatmap.data_diagonal_v2, aes(x, y, fill = value)) +
theme_light() +
geom_tile(color = "grey95",
lwd = 1,
linetype = 1) +
theme(axis.text.x = element_text(angle = 90),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.border=element_blank(),
panel.grid.major=element_line(color='white')) +
ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
subtitle = paste0("Expressed as % of total patients (n = ", nrow(ATC_lv1)," patients)")) +
geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
# scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
guides(fill = guide_colourbar(barwidth = 0.5,
barheight = 20,
title = "%"))
Fig.4 - ATC groups 2nd level
Heatmaps illustrating co-frequencies among 2nd level ATC drug groups for the most prevalent associations from the 1st level ATC drug group heatmap: (A) ATC A and B, (B) ATC A and N, (C) ATC B and N, and (D) ATC B and J. Each cell represents the relative co-frequency of patients using two drug groups concurrently. For instance, in the heatmap (A), the cell “ATC.B01 x ATC.A02” with a value of 0.333 indicates that 33.3% of the total patients (n=2688) are concurrently prescribed drugs from 2nd level ATC groups B01 and A02.
heatmap_2groups <- function(index1 = NA,
index2 = NA,
verbose = TRUE) {
#' This function plots a heatmap of 2 ATC groups at 2nd level
#' inputs: the indexes on the 2nd levels codes for each group
#' ATC.A: 362:377 (16 codes)
#' ATC.B: 391:395 (5 codes)
#' ATC.J: 425:430 (6 codes)
#' ATC.N: 441:447 (7 codes)
ATC_group1 <- BD_chusj[, index1]
ATC_group2 <- BD_chusj[, index2]
# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA,
ncol = ncol(ATC_group1),
nrow = ncol(ATC_group2)))
colnames(heatmap.matrix) <- names(ATC_group1)
rownames(heatmap.matrix) <- names(ATC_group2)
ATC_groups12 <- cbind(ATC_group1, ATC_group2)
for (g1 in seq_along(ATC_group1)) {
for (g2 in seq_along(ATC_group2)) {
heatmap.matrix[g2, g1] <-
nrow(ATC_groups12[, c(g1, g2 +
ncol(ATC_group1))][ATC_groups12[,g1] == 1 &
ATC_groups12[,g2 + ncol(ATC_group1)] == 1,])
}
}
heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value")
# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(BD_chusj), digits = 3)
family.group1 <- substring(names(ATC_group1)[1], 5, 5)
family.group2 <- substring(names(ATC_group2)[1], 5, 5)
p <-
ggplot(heatmap.data, aes(x, y, fill = value)) +
geom_tile(color = "grey95",
lwd = 1,
linetype = 1) +
theme(axis.text.x = element_text(angle = 0),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
ggtitle(label = paste("Relative co-frequences in 2nd-level ATC drug groups",
family.group1, "and", family.group2),
subtitle = "Expressed as % of total patients (n = 2688)") +
geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
# scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
# coord_fixed() +
guides(fill = guide_colourbar(barwidth = 0.5,
barheight = 7,
title = "%"))
# Find most prevalent combinations
if (verbose) {
heatmap.data %>%
arrange(desc(value)) %>%
mutate(perc = round(value / nrow(BD_chusj), 3)) %>%
relocate(y) %>%
rename(!!family.group1 := y, !!family.group2 := x, freq = value) %>%
head() %>% print()
}
return(p)
}# Plot all heatmaps together
p1 <- heatmap_2groups(index1 = 362:377, index2 = 391:395, verbose = F) # A.B
p2 <- heatmap_2groups(index1 = 391:395, index2 = 441:447, verbose = F) # B.N
p3 <- heatmap_2groups(index1 = 391:395, index2 = 425:430, verbose = F) # B.J
p4 <- heatmap_2groups(index1 = 362:377, index2 = 441:447, verbose = F) # A.N
cowplot::plot_grid(p1, p4, p2, p3,
ncol = 2,
nrow = 2,
rel_heights = c(2,1),
scale = 0.9,
labels = c("(A)", "(B)", "(C)", "(D)"),
label_size = 14
)
Fig.5 - ICU-admittance & Deceased
Heatmaps illustrating co-frequencies of exposure among 1st level ATC drug groups for (A) ICU-admitted patients and (B) deceased patients. Each cell represents the relative co-frequency of patients using two drug groups concurrently. Taking the heatmap (A) as an example, the cell “ATC.B x ATC.A”, with a value of 0.934, indicates that 93.4% of ICU-admitted patients (n=858) are simultaneously prescribed drugs from ATC groups B and A. The diagonal cells, such as the cell “ATC.B x ATC.B” showing 0.997, can be interpreted as the relative frequency of the use of that specific ATC drug group alone, signifying that 99.7% of ICU-admitted patients are prescribed only drugs from ATC group B. It’s important to note that the heatmap is mirrored, meaning data below the diagonal mirror those above it. Therefore, the cell “ATC.A x ATC.B” has the same value as the cell “ATC.B x ATC.A”.
# Fig.5A ------------------------------------------------------------------
ATC_lv1 <- BD_chusj[, c(469, 378:390)]
ATC_lv1 <- cbind(IntensiveCare = BD_chusj$IntensiveCare,
ATC_lv1)
# Filter only those patients in Intensive Care
ATC_lv1 <- ATC_lv1 %>% filter(IntensiveCare == 1)
ATC_lv1$IntensiveCare <- NULL
# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA, ncol = length(ATC_lv1), nrow = length(ATC_lv1)))
colnames(heatmap.matrix) <- rownames(heatmap.matrix) <- names(ATC_lv1)
for (i in 1:length(ATC_lv1)) {
for (j in i:length(ATC_lv1)) {
heatmap.matrix[i,j] <- heatmap.matrix[j,i] <- nrow(ATC_lv1[,c(i,j)][ATC_lv1[,i] == 1 & ATC_lv1[,j] == 1,])
}
}; rm(i, j)
heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value")
# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(ATC_lv1), digits = 3)
# store graphic
fig.5A <-
ggplot(heatmap.data, aes(x, y, fill = value)) +
geom_tile(color = "grey95",
lwd = 1,
linetype = 1) +
theme(axis.text.x = element_text(angle = 90),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
subtitle = paste0("Expressed as % of patients in Intensive Care (n = ", nrow(ATC_lv1), " patients)")) +
geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
scale_fill_gradient(low = "#F9FFFD", high = "#a66e29") + #8068b0 #a66e29
guides(fill = guide_colourbar(barwidth = 0.5,
barheight = 20,
title = "%"))
# Fig.5B ------------------------------------------------------------------
ATC_lv1 <- BD_chusj[, c(469, 378:390)]
ATC_lv1 <- cbind(Outcome = BD_chusj$Outcome,
ATC_lv1)
# Filter only those patients in Intensive Care
ATC_lv1 <- ATC_lv1 %>% filter(Outcome == 1)
ATC_lv1$Outcome <- NULL
# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA, ncol = length(ATC_lv1), nrow = length(ATC_lv1)))
colnames(heatmap.matrix) <- rownames(heatmap.matrix) <- names(ATC_lv1)
for (i in 1:length(ATC_lv1)) {
for (j in i:length(ATC_lv1)) {
heatmap.matrix[i,j] <- heatmap.matrix[j,i] <- nrow(ATC_lv1[,c(i,j)][ATC_lv1[,i] == 1 & ATC_lv1[,j] == 1,])
}
}; rm(i, j)
heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value")
# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(ATC_lv1), digits = 3)
# store graphic
fig.5B <-
ggplot(heatmap.data, aes(x, y, fill = value)) +
geom_tile(color = "grey95",
lwd = 1,
linetype = 1) +
theme(axis.text.x = element_text(angle = 90),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
subtitle = paste0("Expressed as % of deceased patients (n = ", nrow(ATC_lv1), " patients)")) +
geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
scale_fill_gradient(low = "#F9FFFD", high = "#cc3746") +
guides(fill = guide_colourbar(barwidth = 0.5,
barheight = 20,
title = "%"))
# Plot all heatmaps together
cowplot::plot_grid(fig.5A, fig.5B,
ncol = 2,
nrow = 1,
# rel_heights = c(2,1),
scale = 0.9,
labels = c("(A)", "(B)"),
label_size = 14
)